Expression Templates for Truncated Power Series

Truncated Power Series technique (Differential Algebra or DA) is a powerful tool for non-linear map analysis of accelerators. The most natural language for numerical DA’s is C++, since it is object oriented and has operator overloading. Traditional C++, though, can be inefficient for scientific programming due to creation of many temporaries and extra loops in overloaded operators. Recent Expression Templates technique allows a user to combine the elegance of object oriented approach with the speed of procedural languages. The way it was created, it is not directly applicable for DA. We created a set of classes whose structure will be suitable for implementing DA vectors and maps. Classes realizing the Expression Templates technique are separated from the client classes, which allows their reuse for different mathematical concepts. Speed tests on KCC compiler showed that new C++ classes for DA have the same speed as hand-coded C. 1 PROBLEMS OF NUMERICAL IMPLEMENTATION OF DIFFERENTIAL ALGEBRA Numerical differential algebra methods (see Ref. [1]) have found increasing use as they can be applied to arbitrary dynamical systems, for example accelerators or other beam and optics devices. The effect of passing of particles or rays through the system can be described mathematically by a map relating the final coordinates (coordinates in the most general sense: momentum, position, charge, mass etc.) to the initial. The map contains the information needed to evaluate such quantities of interest as nonlinear oscillation frequencies and chromatic aberrations. However, except for several trivial cases, it is impossible to find a closed mathematical solution for the map. But one can find solutions up to some order in distance from the central trajectory through perturbation theory based on power series expansion of the map. Numerical differential algebra methods permit one to carry out this perturbation theory for very high order maps, e.g., twelfth order in six variables, where there are millions of separate terms in the map. For the purpose of this paper, one need to know only a couple of facts about DA vectors. In principle, it is an array of coefficients, obtained from Taylor expansion of functions in multidimensional space. These vectors can be added, subtracted, multiplied etc., and the result of most operations depends on the whole set of coefficients of operands. Thus, DA vectors can be represented by a vector-type class with overloaded operators. That is the reason why the implementation of DA methods is easiest within object oriented programming (see [2]), which provides a powerful triad: encapsulation, inheritance, and dynamic binding. C++ is especially well suited since it provides operator overloading. Unfortunately, a big performance price must be paid to realize many of these benefits. Performance losses of a factor up to ten times (exact number depends on vector length) have been reported for codes rewritten in C++. The performance loss arises because usual C++ programming practices lead to the creation of many temporaries and the separation of complicated loops into loops with fewer operations. The latter is especially bad in current times, where RISC chips may have several arithmetic units that can be working simultaneously. As a concrete example, we consider the evaluation of the addition among these vector objects in the following code fragment. Vector y, a, b, c; y = a + b + c; With conventional Vector objects, usually defined as one-dimensional arrays, perhaps including array limits, with properly overloaded arithmetic operators, the typical C++ compiler generates code equivalent to: int n = a.size(); Vector temp1 = new double[n]; for(int i=0; i DA& operator= (const DAExpr& result){ for (int i=0;i static double Add (const A& a, const B& b, int i) {return a[i]+b[i];} }; This class has typical (for vector-like classes) members: double* data and int length, describing the set of coefficients and their number. It has a constructor from int, a copy constructor and a typical destructor. In order to insure a normal behavior of assignment, we need to provide assignment to DA with deep copying (memory management with garbage collection can change its implementation). To implement inlining of all expressions containing vectors, we need assignment to DA expressions. Class for DA expressions is template and relies on behavior of operator[] (int), which is governed by the type of the template parameter A: template class DAExpr { private: A iter; public: DAExpr(const A& a) : iter(a){} double operator[](int i) const {return iter[i];} }; Assignment to DAExpr is a template member function of a no-template class. Not so many compilers provide this capability: KCC compiler (Ref. [4]) is the only one we know. If the compiler does not have template members, see Ref. [3] for the trick Consider how all this should work in a simple example: DA y, a, b; y = a + b; Since this operation will call for assignment of y to an expression (a+b), summation of two DA vectors should return a DA expression, whose type is determined by the performed operation. In this case it is a binary operation of + type (DApAdd) between two double arrays: DAExpr > > operator+(const DA& a, const DA& b) {typedef DABinOp < double*,double*,DApAdd > ExprT; return DAExpr (ExprT(a.data,b.data));} When the resulting expression is assigned to a vector, operator[](int), sitting in the assignment operator, applies this operator to DAOpBin (since it is the filling of DAExpr in this case). Inside DAOpBin::operator[](int), function Op::Apply is called (with Op in our case being DApAdd) template class DABinOp { private: A iter1; B iter2; public: DABinOp(const A& a, const B& b) : iter1(a), iter2(b) {}